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We present detailed numerical simulations and analytical approximations of the propagation 
qq". of nucleons above 10 19 eV in the Local Supercluster, assuming that the ambient magnetic field is 

' turbulent, and its strength 0.01 /iG < B rms < 1 /iG. In such strong magnetic fields, protons in 

| the low energy part of the spectrum, 10 19 eV < E < Ec diffuse, while the higher energy particles, 

with E > Ec propagate along nearly straight lines. The magnitude of the transition energy Ec 
depends mainly on the strength of the magnetic field, the coherence length, and the distance to 
the source; for B rms ~ 0.1 /iG, a largest eddy of length ~ lOMpc, and a distance to the source 
• ~ lOMpc, Ec — 100 EeV. Our numerical treatment substantially improves on previous analytical 

approximations, as it allows to treat carefully the transition between the two propagation regimes, as 
well as the effects due to inhomogeneities expected on scales of a few Mpc. We show that a turbulent 
magnetic field B Ilns ~ 0.1 /j,G, close to equipartition, would allow to reproduce exactly the observed 
spectrum of ultra high energy cosmic rays, up to the highest energy observed, for a distance to the 
source d < 10 Mpc, for the geometry of the Local Supercluster, i.e. a sheet of thickness ~ 10 Mpc. 
' Diffusion, in this case, allows to reproduce the high flux beyond the Greisen Zatsepin Kuzmin cut-off, 

QQ , with a soft injection spectrum j(E) oc E~ 2A . Moreover, the large deflection angles at the highest 

energies observed, typically ~ 10° for the above values, would explain why no close-by astrophysical 
\& counterpart could be associated with these events. 

PACS numbers: 98.70.Sa, 98.62. En 
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The structure and strength of extragalactic magnetic fields are very poorly known at present. These fields likely 
consist of two components: A large-scale component that permeates the whole universe and constitutes a primordial 
• • ■ relic of physical processes in the early universe, and a more localized component which is associated with the structure 
of galaxies, galaxy clusters, and superclusters. For the primordial field only upper limits are known, most notably 
from the limits on the Faraday rotation of radio waves from distant powerful radio sources [DM]. For a field of r.m.s. 
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I. INTRODUCTION 



strength B lms and coherence length A c , this constraint is B lms X^ 2 < 10 _9 GMpc 1 ^ 2 . For the component of the field 
that is uniform on the Hubble scale, this constraint reads B c < 1CP 11 G. Theoretical scenarios for the origin of this 
primordial component, e.g., as relic fields produced by the electric currents present during the electroweak or the 
QCD phase transition ||, or by the Biermann battery effect ||, predict much smaller field strengths, typically of the 
order of 10 -20 G. In contrast, positive detections from measurements of Faraday rotation and synchrotron emission 
exist for the localized component Estimates range from 10~ 9 G up to 10~ 7 G over scales of Mpc, and in some 

cases up to 10 -6 G, e.g., in the core of galaxy clusters M. Notably, a broad ~ 50° region centered on the Local Virgo 
Supercluster seems to be endowed with a magnetic field with a uniform component of strength ~ 1.5 x 10 -6 G [pj, 
and a random component of strength ~ 10~ 6 G |J. The coherence length of random components is usually taken to 
be a cell size of ~ 20 kpc . 

The origin of this component is unclear, the most common explanations being the action of a dynamo that ex- 
ponentially amplifies seed fields of the order of 10 -20 G in the flows of the large scale galaxy structures, adiabatic 
compression of strong primordial fields in the range 1CP 12 — 10~ 9 G, or ejection from radiogalaxies Q or normal 
galaxies ||. 

After several generations of experiments, most notably the Haverah Park p!p| , the Akeno Giant Air Shower Array 
(AG AS A) JL1 12 1 , and the Fly's Eye |p| experiments, the origin and the nature of ultra high energy cosmic rays (UHE 



CRs), with energies E > 10 EeV (lEeV — 10 18 eV) constitute an unsolved problem as well. Nucleons above a few 
tens of EeV are limited in their range to at most ~ 50 Mpc, due to the production of pions on the cosmic microwave 
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background (CMB), known as the Greisen-Zatsepin-Kuzmin (GZK) effect [ fL4| . No compelling astrophysical candidate 
for the source of the highest energy events could be found within this distance scale p5H17|] . Other primaries such 
as nuclei, 7— rays, or neutrinos, pose similar problems jljfl . Theoretical scenarios for their origin include powerful 
radio galaxies [^9[ , cosmological 7— ray bursts 20| , and topological defects left over from phase transitions in the early 
universe that took place at temperatures near Grand Unification Scales 

Magnetic fields above ~ 10~ 12 G on Mpc scales significantly influence the propagation of charged cosmic rays up 
to the highest energies ever observed, namely several 100 EeV, by deflecting and delaying them relative to straight 
line propagation with the speed of light. The absence of a significant correlation of their arrival directions with the 
Galactic plane |2^] and the fact that protons at these energies cannot be confined within the Galactic magnetic field, 
suggest that UHE CRs above ~ 1 EeV most likely have their origin in extragalactic sources. 

We recently developed a Monte Carlo code for the propagation of nucleons in extragalactic magnetic fields. We 
demonstrated that a combination of this code with a maximum likelihood approach allows to extract from observed 
distributions of arrival directions, times, and energies information on the nature of the source mechanism as well as 
on the strength and structure of large-scale cosmic magnetic fields, without taking recourse to any specific scenario 
for the origin of UHE CR (|JJ2^]. This was discussed in a qualitative way in Refs. (|^J2^] and applied to the pairs 
of events recently reported by the AGASA experiment Jl^] in Ref. 28 1 . Our studies were motivated by the fact that 
future large scale experiments |2jJ, such as the High Resolution Fly's Eye |30f| , the Telescope Array |f3l| , and most 
notably the Pierre Auger project should allow to detect clusters of > 20, and possibly more, UHE CRs per 
source, if the clustering suggested by the AGASA results is a real effect caused by discrete sources. The recently 
proposed satellite observatory concept for an Orbital Wide-angle Light collector (OWL) J33| might even allow to 
detect clusters of hundreds of events by watching the Earth's atmosphere from space. 

In our studies we adopted a statistical description of the extra-galactic magnetic field as a Gaussian random field in 
Fourier space with zero mean and a power spectrum with (B 2 (k)) oc k nB for k < 2n/l c and (S 2 (fc)) = otherwise, 
where l c is the numerical cut-off scale. Our studies were thus performed for a "homogeneous" magnetic field, and, 
furthermore, assuming near straight line propagation, i.e. B lms < 10~ 9 G. The subject of the present work is 
to generalize this to the case of more complex structures of cosmic magnetic fields, most notably the possibility 
of a locally enhanced magnetic field associated with the Local Supercluster. Indeed, Ryu, Kang & Biermann |34| ] 
have recently argued that the upper limit on the strength of a homogeneous magnetic field, obtained from Faraday 
rotation measures of distant quasars, becomes B < 1//G in filaments and walls of large scale structure, under the 
very reasonable assumption that the magnetic field scales with the matter density. On the other hand, if the two 
paired events observed by AGASA close to the Supergalactic plane are due to discrete sources and not a statistical 
fluctuation, the field strengths in this structure must be smaller than ~ 10~ 9 G. However, Biermann et al. |3if | 
have argued that due to the accretion flow with embedded magnetic fields towards the Supergalactic plane, weak 
confinement is possible, in which case the AGASA pairs could have originated in unrelated sources. The existence of 
strong magnetic fields in the local Supercluster is therefore an unsettled issue. 

Particles above ~ 10 EeV, that propagate over cosmological distances in magnetic fields as strong as B > few x 
10~ 8 G, enter the diffusion regime in which they turn around several times before reaching the observer. However, 
above a certain energy Eq, whose value depends on the strength of the magnetic field, on the coherence length, and 
on the distance to the source, charged UHE nucleons do not diffuse, and rather propagate in a near straight line. 
This scenario, proposed by Wdowczyk & Wolfendale jj6|, has several attractive features with respect to the observed 
UHE CR spectrum. The transition between diffusion and near straight line propagation could a priori explain why 
the UHE CR spectrum appears flatter above ~ 100 EeV than below, although this trend is not confirmed due to a 
lack of statistics at high energy. Moreover, such strong magnetic fields would result in large deflection angles, even 
at the highest energies observed, and would thus explain why no astrophysical counterpart could be associated with 
these events. Finally, the work of Ryu, Kang & Biermann |Q seems to support the idea that the magnetic field could 
be much stronger within walls and filaments of large scale structure, than in the voids. The spectra of UHE CRs in 
such structures have to some extent been discussed for the diffusive limit in Ref. |37j 

Recently, the diffusion of UHE CRs in the local supercluster has been reconsidered in Ref. pq] . These authors 
combined the diffusion approximation at low energies with the straight line propagation approximation at high ener- 
gies, assuming continuous energy losses and an infinitely thick supergalactic plane. The main difference to Ref. |3(t | 
lies in a better treatment of energy losses, a more careful treatment of the transition regime between diffusion and 
rectilinear propagation, a more quantitative appraisal of the diffusion coefficient, and a more detailed discussion with 
respect to present data. One goal of the present work is to overcome the approximations made in Refs. (3^J|^], i.e, to 
provide more accurate flux calculations at "intermediate" energies where both the diffusion and the small deflection 
approximations break down, and to examine the effect of an inhomogeneous distribution of densities and magnetic 
fields in the local Supercluster. This latter effect is especially important as the distance to the nearest cluster (Virgo, 
d ~ 20Mpc) is already larger than the typical scale of these inhomogeneities for which a few Mpc is a reasonable 
assumption. 
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We describe our approach in Section 2, present results of the numerical simulations in Section 3, and conclude in 
Section 4. Natural units are used throughout, i.e., c = h = 1. 



II. NUMERICAL APPROACH 



A detailed account of our approach was given in Refs. |28|,|24j]. Here, we restrict ourselves to a brief description of 
the most important aspects and the generalizations relevant for the present study. 

A specific scenario is characterized by the magnetic field realization and the position of source and observer. For 
each such scenario, we computed a few 10 3 trajectories of nucleons emitted from a given source in random directions 
within a cone with opening angle of typically 90° and with a flat injection spectrum. This was done by solving the 
equations of motion in the presence of the magnetic fields, taking into account pion production losses by sampling 
the relevant differential cross sections, as well as losses from pair production by protons on the CMB as a continuous 
energy loss process. If the nucleon crossed the sphere of radius d around the source on a cap with an opening angle of 
typically a few degrees, centered on the observer, the nucleon was considered as arriving to the observer. In this case 
its injection energy, and the energy, time, and propagation direction upon arrival, the latter defined relative to the line 
of sight to the source, were recorded. The trajectory of the nucleon does not stop here, however, as, in the diffusion 
regime, one needs to calculate the number density of particles at a given point. Therefore, in our calculations, a given 
nucleon propagates until either the propagation time exceeded 10 10 yr, or the distance exceeded twice the distance 
d between the source and the observer, and, on its way, it can give multiple hits on the "detector" . Since particles 
crossing the cap from cither side were considered as arriving, the maximal deflection angle is 180°. The contribution 
of a particular particle to the flux per area is weighted by the inverse of the absolute value of the cosine of the arrival 
direction relative to the line of sight to the source. If the linear dimension of the cap centered at the observer is 
comparable or smaller than the typical gyro radius of the particle and the coherence length of the magnetic field, 
and the emission cone angle is larger than typical deflection angles, this procedure ensures a realistic distribution of 
arrival energies, times, and directions without too excessive use of CPU time. By varying the cone angles within these 
limits, we verified the independence of the results of the exact cone angle values. Our procedure applies to arbitrary 
magnetic field strengths and deflection angles. 

The resulting distributions of UHE CRs can be folded with the source emission timescale Ts in time and with 
the injection spectrum in energy to obtain distributions for arbitrary Ts and injection spectra. Since we consider 
turbulent magnetic fields, defined by a Kolmogorov spectrum, the spatial configuration of the magnetic field is highly 
anisotropic, as seen from, either the source, or the observer. We thus perform ~ 10 runs for a given set of parameters, 
in order to circumvent, and, at the same, examine, the effect of a particular realization of the field on the spectrum. 

We approximate the structure of the local Supercluster as an infinitely extended two-dimensional sheet of half- 
thickness d s in the x — y plane at a position z s . For the density profile p(z) we assumed a Gaussian of the form 



where c is the density contrast and po is the density at the center of the sheet. It is possible to include a further radial 
dependence of the density profile; however, for the sake of simplicity, we will ignore such dependence in the following 
and take the sheet to be uniform in the x — y plane. 

The random component of the magnetic field is then modeled by using the Fourier transformed Gaussian field 
of strength B Ims with the power law spectrum described earlier, truncated below fc m ; n = 2tt/L, with L the largest 
turbulent eddy, and above & max = 2t:/I c , and multiplying it with the profile f(z). This does not strictly speaking 
obey V • B = 0, but we do not expect this to significantly influence the calculated histograms. For the power law 
index we choose ns = —11/3, corresponding to a turbulent Kolomogov spectrum. 

In addition, the continuity equation and Maxwell's equations allow for a coherent field component B c that is parallel 
to the sheet and whose only dependence is on z and is proportional to f(z). We have verified in our Monte Carlo 
simulations that a coherent component of strength considerably smaller than the random component does not change 
results significantly. For this reason we will neglect a coherent component altogether in the rest of this paper. 

Strictly speaking, due to the high conductivity of the extragalactic medium, the accretion flow onto the sheet causes 
an electric field E = v x B, where v is the velocity field. Because v -C 1, we will, however, neglect this field in 
the present simulations. Other interactions, such as with baryons, are also negligible. We note, however, that solar 
modulation of cosmic rays depends on this electric field, which thus may be necessary to include in future studies. 

In our simulations we usually consider one discrete source located at the origin of the coordinate system, and 
the observer location is characterized by its z coordinate z a , and its distance d to the source; the source itself is 
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characterized by a z coordinate z s . In order to reduce the parameter space, we assume, in all simulations, a sheet 
thickness d s = 5Mpc, a largest turbulent eddy of size L ~ 2d s , and a magnetic field cut-off scale l c = L/10. The 
latter is limited by resolution but does not influence the result as long as l c <C L because the magnetic field power is 
concentrated at the largest scales. For this reason, the effective coherence scale of magnetic fields with a Kolomogorov 
spectrum is given by A c ~ L/(2tt). We further assume power law differential injection spectra Q(E) oc E~ 7 , extending 
up to E = E max = 10 4 EeV, but we note that because of the falling fluxes the results and conclusions are independent 
of maximal energies -Emax for E ma x ^ 10 3 EeV. 



III. RESULTS 



Here, unless otherwise stated, we assume a stationary situation, i.e. such that the emission timescale of the source 
is much longer than the diffusion timescale (typically ~ 10 8 yr). Sources with smaller emission time scales tend to 
produce spectra that are peaked at specific energies, and, therefore, individually provide much poorer fits to the 
combined data above ~ 10 EeV. 



A. Diffusion of UHE CRs in an infinite medium 



We first recall the effect of diffusion in an infinite medium on charged UHE CRs. A similar discussion of this topic 
is given in Ref. [p8[ . Ignoring energy losses for a moment, there will be a transition energy Ec such that particles 
with E ^ Eq will propagate along nearly straight lines, whereas particles with E -C Ec will diffuse. This transition 
occurs where the (energy- dependent) time delay te with respect to rectilinear propagation with the speed of light 
becomes comparable to the distance d. Since in the rectilinear regime, E > Eq, the time delay is given by 

2 d 2 L , / E V 2 / d \ 2 f L \ ( B r ~ 



TB "518^ 6X1 ° UoEevJ U^J UMp^Uo^gJ yr ' (2) 

assuming tlb = —11/3, where r g = E/(eE rms ) is the gyro radius (e is the electron charge), this condition translates 
into 

1/2 / R x / r \ 1/2 



E c ~ 150 — — —— EeV . (3) 

This estimate is confirmed by our numerical simulations within a factor of about 2, for our parameters chosen. For 
E Ec, the observable differential spectrum per solid angle (current density) j(E) is unmodified by the magnetic 
field with respect to the injection spectrum and given by 

The diffusive regime, E < Ec, is usually described by an energy-dependent diffusion coefficient D(E), and the 
expressions corresponding to Eqs.(||) and (Q) read 

T£ "^E)' (5) 

and 

respectively. 

Another way of defining the transition energy Ec is to equate Eqs. ([|) and (^) which gives d = D(Eq). For the 
Bohm diffusion coefficient £> B ohm(E) = r g /3, this results in E c ~ 3.3 x 10 3 (d/10 Mpc)(E rm s/10 _7 G) EeV, about 
a factor 20 higher than Eq. (^) for our parameters. These estimates suggest that the transition regime between 
the analytically tractable limiting cases of diffusion and rectilinear propagation may stretch more than one order of 
magnitude in energy and further emphasizes the need for detailed numerical simulations in that regime. 
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For a diffusion coefficient that scales as D(E) oc E m , one thus obtains j(E) oc E 7 m , in the absence of energy 
losses. From the approximate analytical expression 



3 



/i7r fl(s) dkk2 ( B2 ^y 



where B 2 ma — J °° dk k 2 (i? 2 (fc)}, one expects to ~ 1 for r g > L/(2ir) from Bohm diffusion, and m ~ 1/3 for 
r g < L/(2n) in the case of the Kolmogorov spectrum. As a consequence, the correlation between the time delay te 
and E should switch from, roughly, r oc E^ 1 ! 3 at largest time delays, hence smallest energies, as long as r g < Lf (2tt), 
to te oc E^ 1 in the regime of Bohm diffusion, r g > L/(2ir), and eventually to r cx _E~ 2 at the smallest time delays, 
hence largest energies, in the near straight line propagation regime, E > Eq. As we discuss below, our simulations 
indeed show all three propagation regimes, up to the difference that the diffusing regime with D(E) oc E 1 ' 3 takes 
place when r g < 0.1L/(2n), and D{E) oc E otherwise as long as diffusion holds. 

The above picture is easily derived from analytical calculations. However, such calculations cannot reproduce 
very well the transition between the diffusive regime and the regime in which particles propagate nearly linearly. In 
Figure [l], we show, as the thin solid histogram, the spectrum obtained for UHE CRs propagating in a B lms ~ 5 x 10~ 8 G 
magnetic field, over a distance d — 10 Mpc, in an infinite medium, in the absence of energy losses, and for an injection 
spectrum oc E~ 2 . In the low energy region, lOEeV < E < 70EeV, where Bohm diffusion occurs, one recovers the 
slope predicted by Eq. (g), j(E) oc E~ 2 /D(E) oc £T 3 . At high energy, E > 100 EeV, where particles propagate in 
near straight line, one recovers the injection spectrum, i.e., the spectrum is not modified during propagation. The 
transition between the two regimes indeed occurs smoothly, as would be obtained from a simple interpolation. Our 
simulations allow us to derive the magnitude of the diffusion by equating the ratio of numerical fluxes in the diffusive 
and rectilinear regimes with the ratio of Eqs. (p|) and (Q). It is thus not set a priori; we obtain, for this latter: 
D(E) ~3x 10 32 cm 2 s- 1 (E/l EeV)(B rms /l £tG) , i.e. a factor ~ 10 higher than the Bohm diffusion coefficient. This 
is consistent with the above mentioned fact that we observe Bohm diffusion for r g > 0.1L/ (2ir) and suggests that the 
magnetic field structure in the numerical simulation leads to an effective gyro radius that is about a factor 10 higher 
than the analytical expression r g = E/(eB lms ). 

We now turn to the effects caused by energy losses. Their overall effect in the Monte Carlo simulations is shown 
in Figure [l] as the thick solid histogram, for the same physical parameters as above. The qualitative effect of the 
introduction of pion production and pair production losses, for d c± 10 Mpc, is to change the slope to roughly 
j(E) oc E~ 2 6 (7 = 2), and to deplete the flux at high energies, in favor of the flux at low energies, as pion production 
occurs mainly above the GZK cut-off at ~ 50 EeV. 

Analytically, energy losses are often approximated as being continuous, in which case the rate of change of the 
particle energy is governed by the equation dE/dt — b{E). We note that this is not necessarily a good approximation 
because the dominant energy loss process, namely pion production, is purely stochastic and not continuous. Our 
Monte Carlo simulations do not make this approximation. 

An often used analytical expression for the flux from a point source at distance d in the rectilinear regime |3S|| which 
generalizes Eq. (|J) is then given by (we neglect redshifting) 

Q[Ej(E,d)] dE t {E,d) 
3{E) = ^Wd 2 dE—> (8) 

where Ei(E, r) satisfies the equation d r Ei(E, r) = —b[Ei(E, r)] with Ei(E, d) — E; dEi(E, d)/dE is the ratio of energy 
intervals at emission and observation. 

In the diffusive regime the problem can be described by a partial differential equation for the local density n(E, r) = 
4:TTj(E 7 r) of UHE CRs of the following form 

d t n(E, r) + d E [b(E)n(E, r)] - V [D(E)Vn(E, r)] = q(E, r) , (9) 

where q(E,r) is now the differential injection rate per volume, Q(E) = J d 3 rq(E, r). Eq. (|9) also holds if D(E) is 
space dependent (see next section). In the homogeneous case, an analytical solution of Eq. (§) is known J40| which 
for a point source can be written as 
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E 2 



D(E') 
\b(E')\ 



1/2 



(11) 



denotes an effective path length against losses [X(Ei, E 2 ) = if Ei > E 2 \. 

In the simplest approximation, Eq. (||) and analytical solutions of Eq. @ are just used for E > Ec and E < Ec, 
respectively |56 38 1. For the diffusive solution, the above integral should then be truncated at Eq- A somewhat better 
analytical approximation can actually be obtained by treating particles as propagating rectilinearly as long as E > Eq 
in flight, but diffusively otherwise. For a point source, this amounts to adding in "secondary sources" in Eq. (|^), 
corresponding to the spatial distribution of high energy particles, with E > Eq at injection, that have propagated 
rectilinearly and lost energy down to Eq. The source term is thus replaced by: 



q(E,r) -^q(E,r)e(E c - E) 



ri 



g[£, (£ c , |r - r : | ) , ri ] (£ c , |r - r ; | ) 



47r|r — ri 



dE 



\b(E c )\6(E - E c ) 



(12) 



on the r.h.s. of Eq. (^), where <d(x) is the Heaviside function, n denotes the "secondary" source location; for a 
point-like source, q(E, r) = Q(E)8 3 (r — r s ). This changes the flux towards higher values if Ec is larger than the GZK 
cutoff - 50 EeV. 

We show the analytical predictions Eqs. (||) and (|]) (in the absence of energy losses), and Eqs (|), @, and (§) 
(including energy losses) as light and thick solid lines respectively, in Fig . |l|, for comparison with the Monte Carlo 



results. The continuous energy loss rate b(E) was taken from Refs. 39 41]]. The diffusion coefficient that enters the 



analytical solution of the diffusion equation has been normalized to the value obtained from the numerical solutions 
for a homogeneous medium in the absence of energy losses. We note that the analytical solution, when energy losses 
are included, does not reproduce correctly the transition between the two different propagation regimes. As would 
be naively expected, the energy losses around the GZK cut-off tend to smooth out the spectrum in this transition 
region. Up to this inadequacy of the sudden transition approximation, and the difference in energy losses (obvious at 
high energy), there is good agreement between the two methods, in this standard case. 



B. Diffusion of UHE CRs in the Local Supercluster 



When diffusion occurs in a sheet of finite thickness, the above picture is significantly modified. For a better 
comprehension of the different physical effects that enter the discussion, we distinguish here two different regimes: 
one in which the diffusive regime takes place only for particles with energy less than the GZK threshold ~ 50 EeV, 
and the other in which particles with energies greater than the GZK cut-off also diffuse. For a distance d ~ 10 Mpc, 
and the above choice of parameters, these two regimes roughly correspond respectively to B lms < 5 x 10~ 8 G, and 
B rms >Kr 7 G. 

We first discuss the case where diffusion occurs only below the GZK cut-off. Since the observed UHE CR spectrum 
between ~ 5 EeV up to the GZK cut-off can be fit with one power law, we select B rms ~ 5 x 10 -8 G as a fiducial 
value, in which case particles up to the GZK cut-off diffuse, for a source distance d ~ 10 Mpc. The main effect due 
to the finite thickness of the sheet is that particles can escape freely at the boundary of the sheet, since the magnetic 
field outside is too weak to turn the particles around. Therefore, the energy spectrum of UHE CRs is depleted in the 
energy range where particles diffuse, whenever the source and/or the observer are located near the boundary, or when 
the distance to the source is much larger than the sheet thickness. It is important to note that this effect occurs only 
for those particles that diffuse, as the high energy particles that travel in near-straight line are not affected by the 
presence of the boundary. We recall that diffusion of UHE CRs in the local Supercluster was precisely suggested as 
a mechanism to enhance the local density of low energy particles, with respect to that of high energy particles. We 
thus see that the free escape of UHE CRs at the sheet boundary tends to cancel this enhancement . 

This main effect is shown as the thin solid histogram in Fig. ||, for a source distance d = 10 Mpc, a sheet half- 
thickness d s — 5 Mpc, an injection oc E~ 2 , other physical parameters as above, and no energy losses taken into 
account. The same case with energy losses included is shown as the thick solid histogram. The direct comparison 
of this figure with Fig. [j], in the absence of energy loss, shows that the flux of particles that diffuse in the sheet, is 
depleted by a factor ~ 5, with respect to the flux of particles diffusing in an infinite medium. 

The shape of the spectrum is, however, not significantly changed, and still goes as Q(E)/ D(E). This is in agreement 
with the analytical solution for the sheet geometry, assuming j(E) — at the sheet boundaries. For a point source 
this solution which generalizes Eq. (10) can be written as 
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where r Q , z , r s , z s respectively denote the cylindrical coordinates of the observer and of the source, and g m (r , r s ) is 
given, in the absence of energy losses, by 
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where Kq is a modified Bessel function of order 0, and, when energy losses are included: 
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These analytical solutions are shown in Fig. ^ as light and thick solid lines, respectively. In this figure, the solution 
that takes energy losses into account does not mix the two populations of particles, E > Eq and E < Eq, i.e., 
"secondary sources" are not included here, as they were for the infinite medium. However, as numerical solutions 
of Eqs. @ and ( |l2| ) show, inclusion of secondary sources does not change the result significantly in this case where 
Ec ^5 70 EeV. We note again that in our simulations, the nucleons are propagated until either the distance traveled 
is larger than a Hubble time, or twice the distance from the source to the observer. This means that, in practice, 
the sheet is not infinite, and rather corresponds to a disk, of radius 2d, and half-thickness d s . It is possible to obtain 
analytical solutions corresponding to that case, although such expressions are rather cumbersome [ flO| ]. Moreover, the 
above solution should provide a relatively good approximation, as long as 2d > d s , and we are here in this case. 

Let us comment briefly on Fig. |2|. In this figure, the diffusion coefficient that gives the best fit between the analytical 
solution and the Monte-Carlo is a factor ~ 8 higher than the Bohm diffusion coefficient, to be compared with a factor 
~ 10 in the previous case (Fig. |l|) of the infinite homogeneous medium. This slight difference is due to the fact that, in 
the above analytical solution, the diffusion coefficient is uniform in space, i.e. the solution does not take into account 
the profile of the magnetic field in the sheet. The diffusion coefficient, as averaged over the thickness, is then ~ 1.25 
times larger than the corresponding uniform diffusion coefficient, for the same magnetic field strength in the middle 
plane of the sheet. In the infinite homogeneous case, there is no such profile, and the above factor compensates for 
the difference between the diffusion coefficients. We also note, in Fig. [| that the transition between the diffusive and 
rectilinear regimes does not appear to proceed smoothly. There is a peak around ~ 200 EeV, which showed up for 
most spatial realizations of the magnetic field. We discuss further below this dependence of the energy spectrum on 
the particular realization of the turbulent magnetic field. The comparison of the two cases, with and without energy 
losses, in Fig. ||, shows that this peak is smoothed out by energy losses, as expected. Thus, when the geometry of 
the sheet is taken into account, the above simple analytical solutions do not compare as well to the Monte-Carlo 
simulations, as they did for the homogeneous case. 

We wish to stress the importance of the particular geometry of the sheet, as it precisely counteracts the effect of 
diffusion, in depleting the flux at low energies. Several other effects, which we now discuss, add up to this. 

The profile of the magnetic field in the sheet is such that the field is stronger in the middle plane than on the 
boundaries of the sheet. Therefore, if source and observer are located on opposite sides of the middle plane, diffusive 
low energy particles are here as well depleted with respect to non-diffusive high energy particles. This effect is not 
important as long as source and observer are located within ±d s /2 of the middle-plane, and, outside of this range, 
is dominated by the above effect of free escape do the proximity to the boundary. Therefore, we do not discuss it 
further. 

In the presence of a sheet, the geometry of the problem becomes highly anisotropic. It is all the more anisotropic 
when the distance to the source is smaller or comparable to the largest turbulent eddy, i.e, d < L ~ lOMpc in our 
case, where most of the power of the magnetic field lies. Indeed, particles of different energies follow different paths, 
which, depending on the spatial realization of the magnetic field, may, or may not, drive the particles out of the sheet. 
As far as the UHE CR spectrum is concerned, this effect translates itself in very substantial variations of the flux at 
a given energy, in the diffusing regime, Sj(E)/j(E) ~ 1. It is important to note that this effect does not represent an 
overall variation of the flux, but is differential in energy. In Fig. ||, we show two simulations each consisting of 4000 
particles in total, obtained from the same parameters as above, that correspond to two different spatial realizations 
of the magnetic field. These are shown in solid line. We also show the mean flux, as obtained from the average of 
13 such simulations, i.e. 13 such realizations of the field, in solid histogram. And, we also show the upper and lower 
la deviations from the mean, which show Sj(E)/j(E) ~ 1. The two particular realizations are not extreme in any 
respect, and certainly show that the flux is highly dependent on the particular realization of the magnetic field. This 
strong variation implies that analytical solutions, which solve a diffusion equation that is already spatially averaged 
over the realizations of the magnetic field, are invalidated as long as they claim to solve a one source - one magnetic 
field model. In the limit of high energies E, where the influence of the magnetic field is negligible, the scatter Sj(E) 
should go to zero. The finite scatter around the mean flux visible in Fig. |^ and the subsequent ones is due to the 
finite numerical statistics of 4000 particles per field realization. 
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In the case of an infinite medium, there is also substantial variation of the flux at a given energy; it is not as strong, 
however, as when a sheet is present, for the reasons given above. Our explanation is further supported by the fact 
that, for a simulation where the power index of the magnetic field fluctuations is positive, i.e., where the power is 
concentrated on small scales A c < l c = 1 Mpc, the relative deviation of the flux from the mean is not significant. 

Qualitatively speaking, if one wishes to reproduce the observed UHE CR spectrum with the help of diffusion in the 
Local Supercluster, there are several paths to follow. In either case, one has to find a way to increase the flux of UHE 
CR at the lowest energies < 50EeV. Therefore, one could select a situation in which d < d s , so that the boundary 
of the sheet does not play a significant role. The closest energetic galaxy is Centaurus A, at a distance ~ 5 Mpc. In 
Fig. ||, we show the best fit we were able to obtain for d = 5 Mpc, out of 16 spatial realizations of the magnetic field, 
using an injection spectrum oc E~ 2A in order to match the observations; other physical parameters are as above. All 
other spatial realizations did not give a satisfying agreement with observations, and, even in this case, the fit is only 
barely acceptable. However, we certainly cannot exclude all spatial realizations as consistent solutions, in the present 
case. Our simulations only show that such a solution, if it exists, is a priori unlikely, where the a priori refers to the 
statistical distribution of field realizations. 

Another way to follow would be, loosely speaking, to increase the energy losses so as to increase the flux at low 
energies with respect to that at high energies. One should not increase the distance, however, as free escape would 
then affect the low energy part of the spectrum even more dramatically. This is shown in Fig. [| where the flux is 
shown for a source distance d = 25 Mpc, all other parameters as above. However, if one increases the strength of the 
magnetic field, particles with energies higher than the GZK cut-off start to diffuse. Diffusion increases the distance 
traveled to d + te- Pion production is then affected by the magnetic field as soon as te ~ d, which occurs, for 
E > 50 EeV, for B lms > 10~ 7 G, and d ~ 10 Mpc, as seen from Eq. (|J). Furthermore, a stronger magnetic field results 
in a smaller diffusion coefficient, which, in turn, implies a higher number density of diffusing particles, all things being 
equal. Let us now furnish this discussion with more quantitative arguments. 

In Fig. ^, we show an extreme case, where B Tms ~ 3 x 10~ 7 G, and the injection spectrum oc E~ 2A . There, particles 
up to ~ 200 EeV diffuse, and the effective distance that they travel is ~ 1.4 x d ~ 14 Mpc for E ~ 100 EeV, and 
~ 3 x d for E ~ 70 EeV, with a large relative scatter in the distance traveled of order 1, at all energies. In this 
extreme case, pion production is very strong: particles are strongly depleted at high energies due to energy losses, 
while lower energy particle abundances are increased due to diffusion and energy losses of higher energy particles. 
Because all three regimes in the correlation between te and E described in the previous section appear in this case 
above 10 EeV, we show in Fig. ^, the corresponding distribution of UHE CRs in the time delay-energy plane. There 
one can see the correlation te oc E~ 2 at high energies E > 200 EeV, and the associated scatter, that results from 
stochastic pion production losses. Below these energies, where particles diffuse, the correlation switches to te oc E^ 1 
in the Bohm regime, and to te oc E~ x I 3 for E < 60 EeV. Note that the gyro radius of a 60 EeV proton in a 3 x 10~ 7 G 
magnetic field is r g ~ 200 kpc ~ 0.1L/(2n). Our simulations thus indicate that the D(E) oc E 1 ^ 3 correlation only 
arises when r g < 0.1L/(2tt), and D(E) oc E otherwise. Furthermore, Figs. | and @ indicate that the transition regime 
between diffusion and rectilinear propagation is rather narrow, a factor 2-3 in energy. This is consistent with the two 
estimates of the transition energy Eq from Eq. (||) and from d = D(Ec), if the diffusion coefficient derived from the 
simulations, D(E) ~ lO-DBohm(-E'), is used in the latter. 

The best, and consistent, fit to the data is obtained for the intermediate value of the magnetic field strength, 
Brms — 10~ 7 G, and an injection spectrum oc E~ 2A . Amazingly, this injection spectrum has been derived for UHE 
CRs from clusters of galaxies [Q. In Fig. g, we show the resulting UHE CR spectrum, with the associated upper 
and lower la deviations. Most realizations of the magnetic field give, in this case, acceptable fits to the observed 
spectrum. 

For the sake of illustration, we show in Fig. ^|, the angular images at different energies, as indicated, associated 
with a typical realization corresponding to the case shown in Figs. [2] and |T^, but for our standard injection spectrum 
with 7 = 2.4. Interestingly, a cross-over from several images of the same source at low energies to only one image 
at the highest energies can be observed. This transition occurs at an energy for which the linear deflection d0E (0e 
being the deflection angle) becomes comparable to the field coherence scale A c . This would in turn allow to estimate 
A c . Furthermore, for E > 200 EeV, the typical deflection angle is ~ 13°. This gives the order of magnitude of the 
angular box in which one should look for counterparts associated with the highest energy events observed. 

Finally, for continuously emitting sources, the source power Q required to reproduce the spectra discussed here is 
easily estimated as 



Q = J dEQ(E) = ~ (4n) 2 d 2 E 2 j(E)m^ [{E max / E) {2 ^ , (E min /Ef- 

( I? ^) 2 niax[(i? max /i?) (2 ^ ) ,(i? min /£;) (2 ^ ) ] ergs^ 1 (16) 



. 2 

7.2 x 10 45 ' 
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at energies E where deflection is small. Here, -Emm is the low energy cut-off of the injection spectrum (relevant if 
7 > 2), and the second expression has been obtained from the observed flux at E ~ 300 EeV. This model, for a single 
source, active over > 1 8 yr, thus requires a very high cosmic ray luminosity, in agreement with previous estimates [ p6[ . 
For comparison, the observed luminosity of NGC 4151, a powerful Seyfert galaxy of the local Supercluster, located at 
a distance ~ 13Mpc (for Ho — 75 km/s/Mpc), is L ~ 10 45 erg/s, and the typical bolometric luminosity of a Seyfert 
galaxy is L ~ 10 44 — 10 46 erg/s. On the theoretical side, expectations for the UHE CR luminosity are uncertain and 
vary greatly with the type of accelerator and the emission timescale, e.g., Ltjhecr ~ 10 erg/s for acceleration in 
Active Galactic Nuclei over ~ 10 8 yr, Luhecr ~ 10 51 erg/s over ~ 100 s for fireball models of 7— ray bursts, up to 
^uhecr ^ 10 52 erg/s over > 1 s — 10 4 yr (spread in photons, neutrinos, and nucleons), for topological defects, and up 
to iuHECR ~3x 10 47 erg/s for radiogalaxies Q. 

A one source model, such as considered above, certainly demands a huge output of energy in UHE CRs, when 
one integrates the required luminosity over the required emission timescale. Nevertheless, we note that intermittent 
sources, with pulses of duration T <C te ~ 10 8 yr, with inter-pulse durations AT <C te, and, preferably T <C 
AT, would mimic a continuously emitting source, as the scatter in time delays Ste ~ te ^> AT. Obviously, the 
energy budget would then be greatly alleviated. In fact, the hot spots in radiogalaxies that are believed to be 
UHE CR sources JHJ show intermittent behavior on time scales of ~ 10 4 yr and smaller. Also, nothing prevents us 
from considering more than one source; for instance, there are at least seven powerful Seyfert galaxies in the Local 
Supercluster, whose distance lies between ~ 5 Mpc and ~ 20 Mpc. As well, the combination of many bursting sources, 
with emission timescales <C 10 8 yr, would reproduce a spectrum similar to that of a single continuously emitting 
source. 

Finally, we note that the injection spectrum index that we considered here, 7 = 2.4, is very close to the value 
expected from acceleration in ultra-relativistic shocks ~ 2.2 ±0.1 



IV. CONCLUSIONS 



We have discussed the propagation of charged UHE CR nucleons in a locally enhanced extra-galactic magnetic field 
of strength B lms > 10~ 8 G, in the Local Supercluster. Such high values for the ambient magnetic field in walls and 
filaments of large scale structure obey observational constraints such as Faraday rotation measurements and seem to 
be not unrealistic, as recently advertised by Ryu, Kang & Biermann 34 1. In this regime, UHE CRs diffuse up to 
~ 100 EeV for B rms ~ 10~ 7 G, and a distance to the source d ~ 10 Mpc, and propagate recti-linearly beyond, as first 
advocated by Wdowczyk & Wolfendale |3q ]. As a consequence, for sources emitting on timescales large compared to 
typical propagation times, the spectrum in the diffusive range is steeper than the injection spectrum, and typically 
scales as Q(E)E~ ' 6 , with Q(E) the injection spectrum. This is a result of an interplay between the energy dependent 
diffusion coefficient and energy loss rates. For sources emitting on smaller timescales, the spectrum tends to be peaked 
at a certain energy that depends on the turn-on time of the source. Diffusion models present many advantages in 
the context of UHE CRs. For instance, the switch off of diffusion at high energy might provide a simple way to 
reconcile the seemingly observed flattening of the UHE CR spectrum above the GZK cut-off with a single power law 
injection spectrum that is relatively soft, Q(E) oc E~ 2A . Moreover, the large deflection angles associated with the 
strong magnetic fields could explain why no astrophysical counterpart could be associated with the highest energy 
event, within a ~ 2° error box. 

The problem has been recently reconsidered in Ref. using an analytical solution of the diffusion equation in 
an infinite medium, matched to the near straight line propagation regime at high energy. Here, we have performed 
numerical simulations of the full problem, including a plane geometry model of the structure of the local Supercluster, 
and stochastic pion production energy losses. Our simulations improve on analytical solutions, as they allow to 
simulate the transition between the diffusive and the non-diffusive regimes, and the non-trivial geometry (here modeled 
as a planar sheet), with inhomogeneity length scales for the density and magnetic fields of, typically, a few Mpc. 

Our results differ from the analytical solutions of diffusion in an infinite medium. As a main effect, we find that 
the free escape of particles at the boundary, due to the negligible magnetic field outside of the sheet, depletes the low 
energy diffusive part of the UHE CR spectrum with respect to the high energy (non-diffusive) part. For a distance 
comparable to the thickness of the sheet, this depletion is a factor ~ 5 as compared to the flux predicted for particles 
diffusing in an infinite medium (no boundary). We also find that the highly anisotropic geometry of the problem, 
which is all the more anisotropic when the source distance is smaller or comparable to the largest turbulent eddy, 
implies a strong relative variation in flux at any energy in the diffusing part, of order 1. Furthermore, our simulations 
recover the scaling of the diffusion coefficient indicated by analytical estimates, but with about a factor 10 higher 
normalization both in the infinite medium and for the sheet geometry. 

We thus find that, for sources emitting over timescales larger than ~ 10 8 yr: (i) if the magnetic field is not strong 
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enough, i.e. B lms < 5 x 10~ 8 G, the diffusion coefficient is not small enough to compensate the depletion of low 
energy diffusive particles by ejection out of the sheet; no entirely satisfying fit to the observed UHE CR spectrum can 
be obtained in this case. This result holds even for small distances, down to d ~ 5 Mpc, as the sheet half-thickness 
is only ~ 5 Mpc; (ii) if the source distance is significantly larger than the thickness of the sheet ~ 10 Mpc, or if 
source and/or observer are located near the boundary of the sheet, within ~ ±2 Mpc, there is no consistent fit to the 
observed UHE CR spectrum. Indeed, in this case, pion production becomes significant and depletes the high energy 
part of the spectrum, i.e. above the GZK cut-off, while the diffusive particles feel quite strongly the presence of the 
sheet boundary. 

However, we find a very nice and consistent fit to the observed UHE CR spectrum for a source distance d ~ 10 Mpc, 
and a magnetic field B Ims ~ 10~ 7 G. Smaller magnetic field strengths suffer from the above problem, while stronger 
magnetic field imply too much pion production at high energies. Indeed, for the above magnetic field, particles up 
to ~ 100 EeV diffuse, and diffusion considerably increases the effective distance traveled. We thus obtain a solution 
only in a limited region of parameter space. 

A magnetic field B rms ~ 10 G may not be a coincidence after all, as the equipartition field expected in walls 

of large scale structure, such as the Local Supercluster, is -B cq ui ~ 0.3 pG h wo (T/3 x 10 6 K)^ 2 (p h /0.3p c ) 1 ^ 2 [p],^i|, 
where T and pb denote respectively the temperature and baryon density within the sheet, while p c is the critical 
density, and hi 00 the value of the Hubble constant in units of 100 km/s/Mpc. 

Based on the technique developed here, we will perform detailed studies of spectra and angular images of ultra-high 
energy cosmic ray sources for more realistic source and magnetic field distributions in subsequent work. 

Next generation experiments with their potential to record tens to hundreds of events above 100 EeV will be able 
to detect these signatures and exploit the information associated with them. 
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E [eV] 

FIG. 1. Spectrum of UHE CRs propagating in an infinite medium, in a turbulent Kolmogorov magnetic field of strength 
B rms ~ 5 x I CT 8 G, for a distance to the source d = 10 Mpc, a largest turbulent eddy size ~ 10 Mpc, and a numerical cut-off 
scale l c = L/10. The injection spectrum is Q(E) oc E~ 2 ; the spectrum shown is E 2 j(E), with j(E) the observable differential 
flux, such that the deviation from a horizontal line shows the modification of the injection spectrum. The different cases shown 
are: no energy loss, in thin solid histogram; energy losses, i.e. pair production and stochastic pion production, in thick solid 
histogram (arbitrarily shifted upwards from the no-energy loss case); analytical approximations to each of the two spectra, 
obtained from solutions of the equations discussed in the text, in light solid lines (the upper solid line corresponds to the case 
where energy losses are included, as for the Monte-Carlo results). All fluxes have been normalized to the same total number of 
'detected' particles. Analytical and numerical simulations, with energy losses included, differ at high energy E > 300 EeV due 
to the difference between continuous and stochastic energy losses. 
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FIG. 2. Same as Fig. [I], except that UHE CRs now propagate in a sheet of finite half-thickness 5Mpc; both source and 
observer are located close to the middle-plane. 
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FIG. 3. Spectrum of UHE CRs propagating in a sheet, with the same parameters as in Fig. g The solid histogram shows 
the mean flux, as obtained from the average of 13 fluxes corresponding to different spatial realizations of the magnetic field; 
the dashed lines show the upper and lower one-sided la deviations from this mean. The light solid lines show two different 
spectra, corresponding to two different realizations of the magnetic field. 
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FIG. 4. We show here the best fit obtained, for a magnetic field B rulB ~ 5 x 10 -8 G, a source distance d = 5 Mpc, and 
an injection spectrum oc E~ 2A ; the data points are obtained from a compilation of the Haverah Park pXj , Fly's Eye Jl3| ], 
and AGASA |ll| experiments. The solid histogram shows the mean flux, obtained from the average of 16 different spatial 
realizations of the turbulent magnetic field, while the light solid line corresponds to one, carefully chosen, best fit, realization 
of the field. 
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FIG. 5. Spectrum of UHE CRs propagating in the Local Supercluster, for a source distance d = 25 Mpc, other physical 
parameters as in Fig. ||. The solid histogram shows the mean flux, as obtained from the average of fluxes resulting in 11 
different spatial realizations of the magnetic field; the dashed lines show the upper and lower one-sided la deviations from this 
mean. 
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FIG. 6. Spectrum of UHE CRs, as in Fig. ||, but for a source distance d = lOMpc, and a magnetic field firms — 3 x 10 -7 G, 
averaged over 10 field configurations. 
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FIG. 7. The distribution of time delays te and energies E for the case shown in Fig. ^| The inter-contour interval is 0.25 
in the logarithm to base 10 of the distribution per logarithmic energy and time interval. The three regimes discussed in the 
text, te cx E~ 2 in the rectilinear regime E > 200 EeV, te cx E^ 1 in the Bohm diffusion regime 60EeV < E < 200 EeV, and 
te cx E^ 1 ^ 3 for E < 60 EeV are clearly visible. 
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FIG. 9. Angular image of a point-like source, corresponding to one particular realization of the scenario shown in Figs. ^|and 
^| as seen by a detector of ~ 1° angular resolution, in different energy ranges, as indicated. An E~ 2A injection spectrum was 
assumed. A transition from several images at lower energies to only one image at the highest energies occurs where the linear 
deflection becomes comparable to the effective field coherence length, as discussed in the text. The inter-contour interval is 0.1 
in the logarithm to base 10 of the integral flux per solid angle. 
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